Path-integral Monte-Carlo simulations for electronic dynamics on molecular chains: 

I. Sequential hopping and super exchange 
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An improved real-time quantum Monte Carlo procedure is presented and applied to describe 
the electronic transfer dynamics along molecular chains. The model consists of discrete electronic 
sites coupled to a thermal environment which is integrated out exactly within the path integral 
formulation. The approach is numerically exact and its results reduce to known analytical findings 
(Marcus theory, golden rule) in proper limits. Special attention is paid to the role of superexchange 
and sequential hopping at lower temperatures in symmetric donor-bridge-acceptor systems. In 
contrast to previous approximate studies, superexchange turns out to play a significant role only for 
extremely high lying bridges where the transfer is basically frozen or for extremely low temperatures 
where for weaker dissipation a description in terms of rate constants is no longer feasible. For bridges 
with increasing length an algebraic decrease of the yield is found for short as well as for longer bridges. 
The approach can be extended to electronic systems with more complicated topologies including 
impurities and in presence of external time dependent forces. 
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I. INTRODUCTION 

Electron transfer (ET) processes are one of the funda- 
mental phenomena in complex molecular systems , the 
most prominent one being the primary step of photosyn- 
thesis @, • In the last decade or so much theoretical 
and experimental efforts have been focused on bridged 
molecular systems where the transfer from a donor to an 
acceptor is mediated by a molecular structure connecting 
them 0, IE IE IE El • Particular attention has been paid 
to DNA complexes with to some extent still controversial 
results |10| . Of foremost importance are donor-bridge- 
acceptor (DBA) systems in the stro ngly rising field of 
molecular electronics 0, 0, 0, 0, Il5| . Based on ad- 
vanced methods of synthetic chemistry, nowadays, bridge 
units with specific chemical and physical properties can 
be produced and linked to donor/acceptor species or con- 
tacted with metal leads and thus integrated into electrical 
circuits. 

Since the pioneering work of Marcus and Jortncr 
IB S El the influence of residual vibronic modes and a 
solvent environment on the ET is known to be extremely 
crucial. In fact, the presence of a dissipative surround- 
ing is the prerequisite for a directed (irreversible) trans- 
port across DBA complexes meaning that bath param- 
eters such as temperature and spectral densities have a 
sensitive impact 0, 0| . Accordingly, basically two dis- 
tinct transport channels have been identified. For high 
lying and short bridges quantum mechanical tunneling 
gives rise to a superexchange mechanism associated with 
a characteristic exponentially decreasing yield with in- 
creasing bridge length. If the bridge is sufficiently long, 
however, sequential transfer processes (hopping) domi- 
nate based on thermal activation. Theoretically the com- 
petition between these channels has been investigated 
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based on Redfield theory and simple system+reservoir 
models in Ref. 3 and more recently in Refs. 0, IE 13 • I n 
these studies the mentioned exponential length depen- 
dence as well as an algebraic decrease of the transfer rate 
in case of a sequential mechanism has been found empir- 
ically. This latter behavior has been also derived explic- 
itly from a more rigorous formulation in Ref. . In most 
of these previous studies symmetric DBA systems were 
considered with degenerate donor/acceptor and degener- 
ate bridge states since they already display the essential 
features and make the analysis particularly transparent. 
Additional simplifying assumptions about the nature of 
the thermal bath coupled to the ET have been invoked to 
obtain explicit results since, as for most dissipative quan- 
tum systems [Tsj . a numerical description is extremely 
demanding. 

Indeed, it has been one of the major challenges in re- 
cent years in physics and chemistry to develop efficient 
algorithms to get access to a detailed analysis of these 
systems. As shown first by Feynman and Vernon |l9| 
and in the 1980s by Caldeira and Leggett [2EI, the path 
integral representation is particularly suited as a starting 
point. Namely, for Gaussian bath fluctuations the envi- 
ronmental degrees of freedom can be eliminated exactly 
leading to the reduced system dynamics. In contrast to 
perturbative approaches as e.g. Redfield theory this for- 
mulation has the merit of being applicable down to zero 
temperature, to adiabatic bath modes, and strong fric- 
tion. A direct numerical evaluation of the correspond- 
ing path integral expression for the density matrix is not 
possible though. The reason for that are the interac- 
tions, non-local in time, among forward and backward 
paths generated by the bath degrees of freedom. Typi- 
cally, the range of these interactions grows with An 
efficient numerical approach that relies on a finite inter- 
action range for sufficiently high temperatures has been 
developed by Makri and co-workers |2lJ and applied to a 
variety of systems meanwhile. 

Presently, the only numerically exact approach is given 
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by path integral Monte Carlo (PIMC) simulations (23. It 
has been successfully applied to the spin-boson system to 
study e.g. the influence of nonequilibrium initial bath dis- 
tributions j the role of vibronic and electronic coher- 
ences |23|,|2J], and electronic transfer rates in large ranges 
of parameter space [2^ . For a DBA system with one in- 
termediate electronic site (three-state system) which can 
be seen as simple model for the primary charge separa- 
tion of photosynthesis PIMC results could describe most 
of the experimental key features [2(|. One fundamental 
problem of the PIMC, though, is intimately connected 
with the constitutive role of interferences in quantum me- 
chanical systems. Known as the dynamical sign problem 
|22j | , it turns out that the number of sample paths needed 
to achieve a sufficient signal-to-noise ratio increases ex- 
ponentially with the simulated system time. To partially 
circumvent this drawback various techniques have been 
introduced and shown to stabilize the simulations con- 
siderably m 00 si . 



The goal of this paper is twofold: On the one hand 
the challenge is to present a PIMC approach applicable 
to larger discrete systems and longer simulation times in 
order to lay the basis to treat bridge systems with more 
electronic sites, a complicated topology including impuri- 
ties, and/or in presence of external periodic driving. On 
the other hand, motivated by the above mentioned works 
on symmetric bridges, our intention is to investigate the 
role of the superexchange and the sequential mechanisms 
in ET within a complete system+reservoir model and 
with a numerically exact approach. This way, we avoid 
simplifying assumptions about the bath, the system-bath 
coupling, or the structure of the bath along the molecular 
bridge. Accordingly, our results can directly be compared 
with analytical rate expressions derived within Marcus 
theory golden rule approaches |32j |. or cluster ex- 
pansion methods Surprisingly, such an analysis of 
numerical results within the context of the available the- 
oretical findings is still missing. 



The paper is organized as follows: We start in Sec. [D] 
by defining the system and the bath and introducing its 
path integral representation. In Sec. lllll the connection to 
the description in terms of Master equations is revealed 
and known results for the electronic transfer rates are dis- 
cussed. Then, in Sec. HVI the PIMC used and developed 
in this paper is presented, which leads to a substantial 
improvement in efficiency compared to previous PIMC 
schemes. Only this allows us to study in Sec. [V] first 
the three-state system and then complexes with longer 
bridges in detail, particularly with very good signal to 
noise ratio and larger times. A Conclusion contains our 
main results and some remarks about subsequent work. 
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FIG. 1: Molecular chain for S = 2 1 / , 2 {d = 6) and constant 
nearest neighbor coupling A. 



II. THEORY 

A. Dynamics of the dissipative d-level system 

We investigate the dynamics of a single electron mov- 
ing on a chain of d — 2S + 1 discrete sites, labeled by 
a discrete variable —S<s<S with spacing 1, sep- 
arated by equal distances a but arbitrary energies he s 
(see Fig. QJ. Electronic motion is facilitated through 
tunneling between sites s and s', with a real tunneling 
amplitude A SiS '. The electronic coordinate can then be 
expressed as 



q(t) = a ■ s(t) , 



(1) 



where —S < s(t) < S. The position operator thus is 
equivalent to the spin S operator 



iS,|s) = 



(2) 



with \s) denoting the (orthonormal) localized electronic 
states. In terms of electron transfer, |— S) and IS") repre- 
sent the donor and acceptor, respectively, while the other 
states are referred to as the bridge states. 

The free d-level system (dLS) Hamiltonian can then be 
written as 



H dhS = ?iE z - ns x , 



(3) 



where E 2 describes the energetic distribution of the elec- 
tronic sites according to 



E z \s) =e s \s) 
while S x governs the tunneling, 

s 



E 



|*X*'I + I*'X«I) 



(4) 



(5) 



Of particular importance is the case of nearest neighbor 
tunneling only with 

S x \s) = A._i, a |s-l) + A. jB+ i|«+l> 
for —S < s < S , 
S X \~S) = A_ s ^ s+1 \-S + l) , 
S X \S) = As,s+i\S-l) . (6) 

We already note here that the propagation of the free 
d-level system is most conveniently performed not in the 
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site representation used above, but rather in the eigen- 
state representation of H dL s (for details see Sec. I1V|I . 
This is completely equivalent to the transformation be- 
tween (localized) Wannier states and delocalized molec- 
ular orbitals (Bloch states) in infinite tight binding lat- 
tices. The site representation is particularly suited, how- 
ever, to include the interaction with solvent and vibronic 
degrees of freedom which later on will be integrated out. 
For this purpose, we embed the dLS in a system-plus- 
reservoir model, leading to the total Hamiltonian [18| 



H 



H. 



dLS 



Hr + Hf 



H dhs - aS z J2 c Q X Q + (aS z ) 2 £ — 



E 



(7) 



The residual degrees of freedom are thus modeled by an 
infinite collection of harmonic oscillators, Hb, which bi- 
linearly couple to the position of the electron (Hi). The 
so-called counterterm quadratic in the dLS position op- 
erator aS 2 , which reduces to a physically irrelevant con- 
stant for 5* = 1/2, prevents a renormalization of the elec- 
tronic energy levels by the oscillator bath 0. As dis- 
cussed in detail in Refs. 0, 0, 0, 113 this provides a 
reasonably accurate description of reality for the great 
majority of ET systems. It turns out that for the elec- 
tronic dynamics the environmental parameters enter via 
the spectral density 
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E 



-5(10 — LO a ) 



(8) 



which effectively becomes a continuous function of ui for 
a condensed-phase environment. The Gaussian statistics 
of the environment is determined by the complex- valued 
bath autocorrelation function which for real time t reads 



L(t) = 



\( fe c « x «(*)) (X> Q A Q (0)^) 



i r 

Jo 



smh(npw/2) 



where (3 = l/k B T. For S = 1/2, Eq. Q represents 
the celebrated spin-boson (SB) model, with numerous 
applications in solid-state physics 01 ■ I n the special case 
of a constant nearest coupling, i.e. A S(S / = 5 s '- Si iA for all 
s, s', the described compound reduces to a generalization 
of the SB model, as it in fact then resembles a "spin- 
S'-boson model" 0|. In the following, we repeatedly 
will take advantage of the striking similarity. To make 
connection to (classical) two-state ET theory, we also 
mention the classical reorganization energy hA c \ of the 
Marcus theory 0, , 



A r , = - 



J(u) 



(10) 



For the dissipative electronic dynamics we focus on the 
time- dependent site populations, 



p Sf ,sM = r &{ 



e iHt / n \s f )(s f \e 



mt / n Wi(0)} , (11) 



which are normalized, ^2 S __ s P Sj , lSi (t) — 1, and where 
Wj(0) specifies the initial state of the total compound. In 
most theoretical and experimental works an initial sep- 
aration of the electron and the environment is assumed 
corresponding to an initial density matrix 



Wi(0) 



i)( Si 



-P{H B -sni£) 



(12) 



with the electron held fixed in state \si) and the bath nor- 
malization Zb = Tr{e _/3ffB } assuring the full system's 
density matrix to be normalized for all times. For a trans- 
fer process across the entire chain one typically prepares 
the electron initially in the donor state, i.e. |s<) = |— S) 
0, Q . Above, [i is the dipole moment of the electron, 
and £ denotes the dynamical polarization of the bath 
which is equilibrated with respect to the initial po- 
sition of the electron. By comparing with Eq. , we see 
that u£ — a^2 a c a X a . As pointed out in Ref. |25j, this 
"standard preparation" often used in ET experiments is 
especially suitable for a theoretical description of thermal 
transfer rates. 

The present formulation can easily be extended in or- 
der to describe the impact of an external time dependent 
driving force. This as well as abandoning the restriction 
to nearest neighbor coupling will be the subject of a sub- 
sequent paper. 



B. Path integral representation 

The path integral representation provides a formally 
exact expression for the dynamics of the electronic pop- 
ulation and is thus the starting point for a numerically 
exact Monte-Carlo (MC) scheme. The standard proce- 
dure is then to eliminate the bath degrees of freedom to 
arrive at the reduced dynamics. This can be done ex- 
actly due to the harmonic nature of the bath. As shown 
in Ref. 0, one thus obtains for Eq. (fTTft 



Ps 



,(*) 



Vs S S {t),s f exp <^ -5 d Ls[s] - $[s_ 



(13) 

Here the path integration runs over closed paths s(t) 
starting at 5(0) — Si and propagating along the real-time 
contour t G — > t — > 0, which connect the forward and 
backward paths s(t') and s'(t'), respectively. Further- 
more, SdLs[s] denotes the total action of the free dLS, 
and the influence of the traced-out bath is encoded in 
the Feynman- Vernon influence functional $[s] |Tflj . In 
terms of the bath autocorrelation function ©, it reads 
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dt' 
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dt"[s(t')-s'(t')\ [L(t' -t")s(t") 
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where 



L*(t' -t")s'(t")} 

\ l dt'[s 2 {t')- s l2 {t')] , 

o 



3 , A") OA 

u= - I duj— — = 2A c i 

7T ,/n UO 



(14) 



(15) 



The influence functional introduces long-ranged nonlo- 
cal interactions among the dLS paths so that in general 
an explicit evaluation of the remaining path integral in 
Eq. I|13fl is possible only numerically. Before we explain 
the details of the PIMC scheme in Sec. IIVI however, we 
focus in the next Section on the description of the popu- 
lation dynamics in terms of Master equations with associ- 
ated transition rates. Corresponding results will provide 
the appropriate tools in order to analyze and understand 
the MC data. 



III. POPULATION DYNAMICS 

A. Master equations 

As has been derived in Ref. [3i| , the dynamics of the 
population for a spin-S'-boson system can be cast exactly 
into the form of a retarded non-local (in time) Master 
equation 



with time dependent transition rates T s ^ Sk (t) the Laplace 
transforms of which are given by appropriate cluster func- 
tions. Provided the dissipative influence is sufficiently 
strong or temperature sufficiently high, the electronic dy- 
namics (|16|) can be simplified to a conventional Master 
equation. In essence, one invokes a coarse graining pro- 
cedure in time so that after some initial transient motion, 
i.e. for t > Ttrans, the population dynamics obeys 



P, 



.(*)= E 



(17) 



--S 



Here, the nearest neighbor transition rates r s . s ±i are 
related to sequential hopping processes, while T-s,s 
describes direct transfer from donor to acceptor with- 
out any real populations in the bridge states (superex- 
change) . 

In the case studied here in more detail, namely, de- 
generate donor and acceptor energies (e_s = £5 = 0), 
degenerate bridge energies (e_s+i = •■■ = f-s-i = e ), 
and nearest neighbor tunneling with A StS i = <5 s '_ S) iA, 
the number of independent rates shrinks considerably. 
Eq. I|17|) becomes particularly simple in case of purely 
sequential transport where it can be written as 



P = AP 



(18) 



dt 



E 



dt' f s , Si (t - t') P s , Si (t') (16) with the rate matrix 
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r_DB (Xbd + r#) 
v B 



V 





-2T B 











-2T B 

r B 




r B 

-(T B + T B a) Tab 
Tba -Fab I 



(19) 



Here, P = (P D , P B i, • • • , Pet, Pa) where P D (P A ) refers 
to the donor (acceptor) population P-s (Ps) and Pbi 
to the population Pi-s of the b — d — 2 bridge states. 
The rate constants Tbb (Pb.d) govern the sequential 
forward (backward) transport between donor and first 
bridge state, Tab (Tad) the sequential forward (back- 
ward) transport between acceptor and last bridge state, 
and Tb the sequential transport on the bridge. Coher- 
ent transport (supcrexchange) would be accounted for 
by additional rate contributions in the upper and lower 
triangular of A. 

According to Eq. i|18|) the populations P s (t) re- 



flect multi-exponential dynamics, with the corresponding 
rates given by the eigenvalues of the rate matrix A. The 
approach to equilibrium happens to be on the relaxation 
time scale tr which is the inverse of the smallest of these 
eigenvalues T^. Forward and backward rates in A are 
connected via detailed balance, 



r 



BD 



pa 
E 

P? 



Tdb — 



2PZ 



1-bPl 



DB 



(20) 



where P^P and P^ denote the equilibrium occupation 
probabilities of the bridge states and the donor, while 
for symmetry reasons PJ 3 = Pjj, Tba = Tbd, and 
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Tab = Tub- The equilibrium occupation probabili- 
ties are Boltzmann distributed with respect to an effec- 
tive electronic parameter e(/3) which, for strong damping 
and/or higher temperatures, converges with the physical 
bias e 34]. It is then straightforward to calculate P s (t) 
from Eq. (|18(l or the more general expression (|17|l which 
in the appropriate parameter range gives an accurate de- 
scription of the exact dynamics for times t larger than 
some transient time scale itrans- The corresponding ex- 
pressions for d = 3 are given in App. [5] 

Hence, once MC data for the populations are gener- 
ated, the associated rate constants are obtained by fitting 
the electronic dynamics obtained from Eqs. (|18|l and l|19() 
to the numerical data. At sufficiently low temperatures, 
however, coherent oscillations govern the population dy- 
namics over a long period of time and the transfer can no 
longer be characterized by time independent decay rates. 
Outside this range information about the impact of su- 
perexchange can be extracted by comparing fits to the 
sequential (|18|l and the complete model l|17|) . 



with the damping kernel Q(t) specified in Eq. i|34|) . The 
golden rule rate expression reduces for high temperatures 
to the nonadiabatic limit of the Marcus formula pijl and 
captures nuclear tunneling at low T. Within a pertur- 
bative treatment in powers of A 2 ™, superexchange rates 
appear as higher order contributions n — 2,3,4,.... In 
case of a three-state system with en = £a, is ob- 

tained in leading order n = 2 as 33] 

r s , GR = 2A 4 Re|^ d Tl dT 2 dT 3 exp[Q*(r 3 ) 

" Q(T 2 + T 3 ) - Q( n + T 2 + T 3 ) - Q(T 2 ) 

- Q(n + r 2 ) + Q(n)] 



x exp[— ie(ri - r 3 )] 



(23) 



Even for this more involved expression simplifications 
arise in certain limits. In the high temperature range 
hf3uj c <C 1, where Q(t) can be expanded for short times, 
one obtains the classical superexchange rate [33| 



B. Sequential and superexchange rate expressions 

As already noted, basically two transport channels 
have been found in previous studies on bridge mediated 
ET, namely, a sequential channel and a superexchange 
channel. The latter one is expected to dominate for high 
lying and sufficiently short bridges, whereas in other cases 
the former one prevails. Here we briefly collect available 
analytical results for the transitions rates, particularly 
for symmetric DBA systems (degenerate donor and ac- 
ceptor energies, degenerate bridge states), to analyze our 
PIMC results accordingly (see Sec. IVj) . 

In the classical limit h/3uj c <C 1 the sequential forward 
rate is given by the well known Marcus result Q 



A 4 



A 2 



DB, Marcus 



Tlfl 



-f3F*(e) 



1 + ttA 2 /(A c1 lu c ) V A cl fc B T 

with the activation free energy barrier F*(e) = h(e + 
A c i) 2 /(4A c i) where e = e# — is the energy gap between 
bridge and donor (acceptor). As long as a rate descrip- 
tion is valid at all, this expression is applicable from the 
weak adiabatic domain (uj c < A), where it eventually 
becomes independent of A, to the nonadiabatic range 
(uj c ^> A). Recently, for adiabatic to moderate nona- 
diabatic environments an extension of this rate expres- 
sion to lower temperatures where quantum fluctuations 
and nuclear tunneling come into play has been derived in 
Ref. [13. 

In the quantum domain Ti(3uj c <; 1 and for small tun- 
neling amplitudes A <C co c (nonadiabatic range) transfer 
rates are calculated perturbatively by invoking Fermi's 
golden rule. For the sequential forward rate the result is 



T DB ,GR(e) = A 2 



dt exp[— iet — Q{t)\ 



(22) 



SCL.GR 



2(e-A cl ) 2 V A c i 



I3F+ 



+ C 



-(3F- 



(24) 



with the activation energies F + = h(e + A c i) 2 /4A c i and 
F_ = h(e - 3A c i) 2 /4A c i. This rate is well-defined even 
in the resonant case e — -> A c j where, however, also the 
incoherent transfer is extremely efficient with a strong 
bridge population. 

In presence of a quantum bath Ti(3lo c •> 1 the classical 
reorganization energy drops out of the golden rule rate, 
see Eq. (|24f) . since nuclear tunneling is relevant. The su- 
perexchange expression can be simplified for a high-lying 
intermediate state where the DBA system can effectively 
be treated as an AD system with effective coupling [33J, 
i.e., 

r S QM,GR~^Y I rfrexp[-4Q(r)] . (25) 

^ J —OO 



IV. SIMULATION METHOD 

As already addressed, the time evolution of a dissipa- 
tive quantum system can in general not be evaluated an- 
alytically. In the past the PIMC method has been proven 
as a very promising approach to obtain numerically ex- 
act results even in regions of parameter space where other 
approximative methods fail. In particular, in recent ex- 
periments on molecular wires contacted with metal junc- 
tions [3r3 ] low temperature measurements down to 30K 
have been performed which necessitates the inclusion of 
strong quantum effects also in the environment. 

For this purpose we turn to the path integral repre- 
sentation (|13|l and employ a discretization described in 
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detail in Refs. |25l l'26| . For the forward and backward 
paths the time axis is sliced via q uniformly spaced points 
with discretization steps r = t/q. The path integral in 
Eq. H13|) then becomes 



F, M (i) = ^^ 1 ,, / /)[{sj}] 



(26) 



with 



P[{*j}] = 



2q 



J| K(s k ,s k+1 ) 



,fc=i 



(27) 



The sum runs over all realizations of the discretized 
spin path {sj} = {si = Sj, s 2 , . . . , s 2g , s 2q +i = sj, and 
K(sj, Sj+i) denotes the coordinate representation of the 
free dLS propagation over the time interval r, i.e., 



K(s,s',t) = (s\enp(-iTH dhS /h)\s') 



(28) 



While in principle for the dLS Hamiltonian this 
propagator can be obtained by exploiting the 
symmetric Trotter splitting exp(— irHdhs/f 1 ) = 
exp(-irE 2 /2) exp^rS^) exp(-irE z /2) + 0(r 3 ) and 
evaluating (sj+i exp(irS x )\sj), a more accurate ap- 
proach follows from the eigenstate representation 



HolsIM = E N \4> N ) , N = 
This way, one obtains 

d 

K(s,s',t) = ^ ( s l<^v}(<^v| s ') c 
N=l 



(30) 



which can be easily computed numerically once the dLS's 
parameters are specified. Another advantage of this for- 
mulation is that it can be immediately applied to all 
quantum systems which can be mapped effectively onto 
a dLS, e.g. by a proper reduction of its Hilbert space. 

To arrive at a discretized form of the influence func- 
tional {SJ, the sum and difference coordinates 

V (t') = s(t') + s'(t') , £(t) ee s(t') - s'(t') (31) 

are introduced, which read -q(t') — rjj (£(£') = £.,-) for 
t' G [(j — 1)t— t/2, (j— 1)t+t/2] in their discretized form. 
The sum paths are also considered as "quasi-classical" , 
while the difference paths capture quantum fluctuations 
[l8L |26]| . Equation l|14f) finally can be written as 



$[*i,»7,£l 



(Si) 



(32) 



i=2 



with 



j>k=2 



= 2 Si I m {Q((i-2)r)-Q((i-l)r)} 



for 2 < j < q , 
A = Rc{Q(r)} , 
X = Im{g(r)} , 

A, = Re{Q((l-l)r) + Q((l + l)T)-2Q(lT)} 

for 1 < I < q - 2 , 
X ; = lm{Q((l-l)T) + Q((l + l)T)-2Q(l T )} 

for 1 < / < q - 2 . (33) 

Here, one has introduced the twice-integrated bath au- 
tocorrelation function Q(t) defined by Q(t) — L(t), 
Q(0) = with <2(0) = iA c j. Explicitly it is gained from 
Eq. 10 to read 

1 f°° J(uj) 
Q(t) = - / dw-V{coth(S, ) Sw/2))[l-cos(wt)] 
~ /o 



+ i sin(u;i)} 



(34) 



In the sequel we consider a spectral density of the form 

J{lo) = 2TiaLoc- u/u - , (35) 

which is equivalent to ohmic damping with a cut-off fre- 
quency u> c . In this case Q(t) can be calculated analyti- 
cally and one obtains 



(29) Q(t) = 2a 



ln(l+iw c i) - In 



r(o + it/hp)Y(tt - it/np) 



(36) 

with Q = 1 + l/(h(3ujc) and the Gamma function T(z). 

Equations l|26|) . (|27|l and l-il'l) constitute a discretized 
form of the populations (fTTl) and thus provide a starting 
point for PIMC simulations. Unfortunately, this method 
is handicapped by the dynamical sign problem |22| . It 
originates from quantum interferences between different 
system paths {sj}, causing a small signal-to-noise ra- 
tio of the stochastic averaging procedure: The exponen- 
tial increase of the paths' configuration space with the 
maximum real time under study results in an exponen- 
tial decrease of the signal-to-noise ratio. Various proce- 
dures to mitigate the sign problem have been proposed 
in the past, like maximum entropy methods |2jfl , filter 
tec hniq ues |28l l29l | or the multilevel blocking approach 
|25l l3l| . Here we employ a filter technique optimized for 
dissipative spin systems as suggested by Egger and Mak 
|26| which exploits the special symmetries of the influence 
functional. 

The backbone of this approach is the observation that 
the absence of a self-energy like term for the quasi- 
classical paths in Eq. (|32|) allows to express the 
summation over the latter as a series of q — 1 matrix 
multiplications. This usually requires drastically less 
computational effort than performing the correspond- 
ing sums in Eq. (|26[) and thus can be carried out ex- 
plicitly. It reduces the degrees of freedom from the 
2q — 1 variables {r]2<j<q+i, £,2<j<q} to the q — 1 quan- 
tum variables {£,2<j<q} and therefore significantly im- 
proves the numerical stability of the corresponding MC 
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simulations. Explicitely, after switching to the quasi- 
classical and quantum coordinates H31|) and defining the 

[2(5'-|^l) + l]x[2(S'-|^+i|) + l] matrices according 
to 

= exp -irjj+i ^2 €k x k-j-i 

\ k>j 

xif(fe+0)/2,fe+i+0+i)/2) 

xK*{{ m - fc)/2, ( Vj+1 - 0+0/2) , (37) 

the population (|26(l becomes 



{&}>f«+l 



fc>j 

xtf^O, (38) 

Here, the sum over the {%< q } has been written as a 
product over the q matrices where the first (last) 

of them is only a row (column) vector due to 771 = 2sj 
(rj q+ i = 2sy). This matrix product can be performed 
explicitely such that only the sum over the quantum 
variables {£,j}2<j<q has to be evaluated by PIMC sim- 
ulations. The resulting mitigation of the sign problem 
allows for numerically stable simulations which consume 
significantly less CPU time than corresponding ones eval- 
uating the standard expression l|2t)|) . Moreover, another 
considerable speedup can be gained by optimizing the 
MC weight with respect to the dissipative regime (see 
App.E). 

We note in passing that the above described ap- 
proach also resembles the multilevel blocking strategy 
[30| : On the first level, the harmonic bath degrees of 
freedom are integrated out, while on the second level the 
quasi-classical coordinates get eliminated; the addend in 
Eq. (|38|l then corresponds to the respective level-2 bonds. 
It furthermore seems noteworthy that, while here we 
only present results for electronic systems with constant 
nearest-neighbor coupling, A SiS / — <5 s '_ Sl iA, the above 
outlined approach is suitable for arbitrary electronic sys- 
tems as long as they can be described by a finite Hilbert 
space. Corresponding applications will be shown in a 
subsequent paper. 



V. NUMERICAL RESULTS 

Next we present results for a symmetric donor-bridge- 
acceptor (DBA) system (degenerate donor and acceptor 
energies, degenerate bridge states with energy gap he be- 
tween bridge and donor/acceptor) with constant nearest 
neighbor coupling obtained from PIMC simulations as 



described above. Albeit its simplicity this model cap- 
tures essential features of bridge mediated ET and can 
serve as a basis for more elaborate studies. In part, this 
work was motivated by former theoretical approaches on 
symmetric DBA systems 0, H, E| which, however, re- 
lied on simpler system+reservoir models and/or approx- 
imate numerical methods. In particular, in contrast to 
these latter, the formulation 0) takes the dynamics of 
the vibronic structure of the DBA compound fully into 
account. Further, the MC approach is not plagued by 
the limitations of Redfield theory, but rather applies also 
to low temperatures, slower bath modes, and stronger 
dissipation. Hence, it reproduces in the proper limits the 
known analytical findings specified in Sec. IIII 51 

Clearly, the ET process considered here eventually 
leads to thermal equilibration of the electronic sites. Due 
to the relation between relaxation rates and conductances 
in stationary non-equilibrium [T^ |. our results give also 
insight, at least qualitatively, into ET through metal- 
molecule-metal junctions. In fact, to get quantitative re- 
sults, the present approach can in principle be extended 
to the latter case by eliminating the electronic states in 
the metal contacts as an additional bath. Corresponding 
work is in progress. 

The superexchange mechanism is a truly quantum me- 
chanical coherence phenomenon. Thus, in order to inves- 
tigate this scenario in more detail within the MC ap- 
proach, we choose parameters such that for lower tem- 
peratures the population dynamics exhibits coherent os- 
cillations, while for slightly higher T, with a bath that 
is still quantum mechanically, a rate description is ap- 
plicable. We take a damping coefficient a = 1/4 well 
below the coherent-incoherent transition a = 1/2 [Isj 
and a moderate cut-off frequency oj c /A = 5 so that adi- 
abatic effects are expected to show up. Note that the 
cutoff frequency oj c also defines the maximum of the spec- 
tral density distribution l|35|l . such that even frequencies 
somewhat larger than co c contribute. The classical reor- 
ganization energy follows as A c i = 2aoj c — 2.5A. Ac- 
cordingly, we expect superexchange, if important, to be 
particularly pronounced. 

To confirm the proper choice of the above parameters 
and to fix the temperature range where a description 
in terms of rate constants applies, we start by present- 
ing results for the equilibrium bridge population of 
a three-state system at a fixed bridge energy e/A = 5 
for varying temperature. As can be seen by comparing 
obtained from imaginary-time MC simulations with 
a simple Boltzmann distribution, i.e. l/[2 exp(/3e)+l], de- 
viations occur for frf]A > 0.3 related to the strong influ- 
ence of quantum derealization (see Fig. 01 . This triggers 
the tendency of coherent oscillations in the population 
dynamics which becomes especially obvious for an ener- 
getically degenerated bridge (see Fig. |3J) . Damped oscil- 
lations are clearly observable in the intermediate state for 
7i/3A > 0.3. These results verify that (i) with the above 
choice of parameters electronic coherent effects are cru- 
cial and (ii) that the ET dynamics can be captured by 
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q\ , 1 , 1 , 1 , 1 , U 

u 0.1 0.2 0.3 0.4 0.5 

npA 

FIG. 2: Equilibrium occupation probabilities of the 
bridge state for S = 1, e/A = 5, a = 1/4 and lu c /A — 5 as 
a function of inverse temperature. Circles denote imaginary- 
time QMC results (with error bars smaller than symbol size) 
while the solid line refers to l/[2e R/3E + l] according to a Boltz- 
mann distribution. 

transfer rates for H/3A < 0.3. With respect to the bath, a 
typical temperature h/3A = 0.1 corresponds at the max- 
imum w max = uj c of the spectral density to ft/3w max = 0.5 
so that a high temperature approximation does not ap- 
ply. The simulations have been performed for system 
times long enough and sufficient stochastical accuracy to 
allow for studying details of the dynamics even close to 
the approach of equilibrium. 

Having fixed the proper range of parameters, we now 
turn to a detailed analysis of the population dynam- 
ics. First, we consider a three-state system with varying 
bridge energy < e/A < 20 and at a fixed tempera- 
ture h(3A = 0.1 (Fig. gj. For bridges with e/A < 7.5 
the intermediate state approaches thermal equilibrium 
within the simulation period. Interestingly, a bridge en- 
ergy e/A = 2.5 leads to an increase of the acceptor popu- 
lation almost as fast as for the degenerate case e/A = 0. 
This is due to the fact that for the former bridge the equi- 
librium population P4 is larger while the bridge is still 
not high enough to considerably reduce the donor-bridge 
rate Tob- The population dynamics for this system can 
now be analyzed with the full i|17|) and the sequential H18|) 
model to extract the rate constants for times t > itrans- 
The rate constants were then obtained by repeating this 
procedure for successively increasing gratis such that the 
corresponding rates eventually saturate at plateau val- 
ues. The latter are depicted in Fig. [S] where also the 
nonadiabatic classical Marcus (|21|l and the golden rule 
(|22|l results are shown. The numerical results basically 
follow the latter ones with minor deviations that can be 
attributed to adiabatic effects (e.g. recrossing phenom- 
ena) not included in the golden rule rate expression. Se- 
quential rates are extracted for all three electronic states 
where those obtained from the donor and acceptor co- 
incide while rates extracted from the bridge exceed the 
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FIG. 3: Donor, bridge, and acceptor populations Pn(t), 
Pfl(i), and PA(t), respectively (from top to bottom), for 
5 = 1, e/A = 0, a = 1/4, uj c /A = 5, and (top to bottom 
for donor, bottom to top for bridge and acceptor) h/3A = 
0.05, 0.1, 0.2, 0.3, 0.4, and 0.5. Error bars correspond to one 
standard deviation in each direction. Dotted lines denote the 
equilibrium population P^ = PJf = Pa =1/3. The solid 
lines are guides for the eye only. 
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FIG. 4: Donor, bridge, and acceptor populations Pn(t), 
Pfl(i), and PA(t), respectively (from top to bottom), for 
S = 1, a = 1/4, uj c /A = 5, h/3A = 0.1, and (top to bot- 
tom for bridge and acceptor, bottom to top for donor) e/A = 
0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, and 20. Error bars correspond 
to one standard deviation in each direction. Dotted lines de- 
note the respective Boltzmann equilibrium populations of the 
bridge state. The solid lines are guides for the eye only. 




FIG. 5: Sequential transfer rates Ydb obtained from the 
bridge state (open diamonds) and donor and acceptor state 
(closed diamonds) for S = 1, a = 1/4, u c /A — 5, and h/3A — 
0.1 as a function of e. The solid line refers to the golden rule 
rate 1221 . the dashed one to the Marcus rate 12 H . 



former by about 5-10%. This deviation signals a short- 
coming of Eqs. (|18fl and 119|) in this parameter regime. 
In fact, the reduction of Eq. I|16|l to Eq. 1|17|) is in a strict 
sense only justified in the strongly damped regime. The 
impact of superexchange, however, is very weak and its 
corresponding rates are smaller by two orders of magni- 
tude. Accordingly, the sequence of superexchange rates 
for various itrans obtained from the fitting procedure de- 
scribed above does not saturate but rather strongly fluc- 
tuates such that a quantitative description is out of range. 

This behavior can be understood from the analyti- 
cal rate expressions. Namely, upon closer inspection 
of Eqs. (|24|) and 121|l one derives that in the classical 
regime Tscl,gr dominates against the sequential trans- 
fer Tub, Marcus whenever hf3e ^> 1 and e » A c i such that 



rscL,GR A 2 



r 



DB, Marcus 



exp(/3^e 2 /4A c i) > 1 



(39) 



This verifies the expected and well-known fact that for 
a classical bath the superexchange mechanism prevails 
whenever thermal activation from D to B is suppressed. 
For a quantum bath, however, the situation is more 
restrictive. From Eq. H23[) one estimates Fsqm,gr °c 
A 4 /(u! c £ 2 ) for h(3oj c > 1 (and a of order 1) such that 
for fixed temperature hf3A < 0.5 superexchange trans- 
port prevails only for extremely high bridge energies, 
i.e. e/A > 40, where the transfer is essentially frozen 
on the time scale accessible in MC simulations and in 
realistic experiments. With electronic couplings of the 
order of 100 cm -1 0,0 , the above condition requires at 
T ~ 150 K bridge energies of the order of e « 5000 cm -1 
leading to typical transfer times of the order of 30 /is. 
Having in mind that our primary focus is to identify re- 
gions where one of the transfer channels, sequential or 
superexchange, dominates, superexchange turns out to 



in 



be negligible in Fig. 

At least as a minor effect we find superexchange in 
the ET dynamics shown in Fig. 0] for e/A > 5. We 
thus performed simulations for fixed bridge energy, but 
varying temperature 0.05 < hf3A < 0.5, see Fig. 
As already mentioned, for hfiA > 0.3 a rate descrip- 
tion is questionable as the bridge population shows signs 
of coherent oscillations. Sequential rates are extracted 
according to the above scheme from donor and accep- 
tor as well (cf. Fig. 0). The results are in qualitative 
agreement with the golden rule description with devi- 
ations between donor and acceptor rates appearing in 
the low temperature range where quantum coherences 
exist. All numerical rates are below the nonadiabtic val- 
ues, again due to adiabatic effects in the bath dynam- 
ics. Superexchange cannot be found. This is under- 
stood from the above discussion: For fixed e the ana- 
lytical superexchange rate exceeds the golden rule one 
only in the temperature range h(3A 3> 2 ln(e/A)/(e/A) 
corresponding to low temperatures hf3A ^> 0.6 where, 
as seen above, for weaker dissipation (a < 1/2) a rate 
description is in a strict sense no longer feasible. The sit- 
uation changes when donor and acceptor are biased with 
an energy gap en A = (d — > 0. Then, based on the 
generalization of Eq. (|23l for cda ^ |3^| one can show 
for a high lying intermediate state that superexchange 
exceeds sequential transport already at higher tempera- 
tures HP A > 21n(e/A)/[(e_ Dj4 + e)/A]. In fact, for biased 
ET systems indications of superexchange have already 
been seen in previous MC simulations for three-state sys- 
tems nu. 

Let us now turn to longer bridges with b = d — 2 > 1 . 
The stability of the MC procedure allows to extract rate 
constants for systems up to b = 10. In Fig. |S] results for 
b = 1 up to b = 10 bridge states are displayed. Appar- 
ently, the donor dynamics basically saturates for b > 2 
on the time scale of the simulations meaning that the 
donor decay depends on the bridge length only at very 
large times and then only slightly. This amounts to the 
fact that the equilibrium populations decrease only al- 
gebraically with b. In contrast, the time evolution of 
the acceptor and its closest bridge states are strongly af- 
fected by the varying bridge length. The numerical data 
can be very precisely reproduced by a sequential trans- 
fer model according to Eq. I|18[) where due to symme- 
try and detailed balance only two fit parameters enter, 
namely the forward rate from the donor Tjjb and the 
bridge rate Tb- It turns out that by fixing these param- 
eters at Fob — 0.284 and Tb = 0.348 all populations 
independent of the bridge length can be described, thus 
verifying that the transport is more or less completely 
sequential. The nonadiabatic golden rule formula 122|) 
predicts a little bit higher values Vdb,gr = 0.297 and 
^b,gr = 0.37 in agreement with our findings for the 
three-state system. Our main focus lies on the behavior 
of the smallest eigenvalue of the rate matrix A, i.e. the 
relaxation rate T^, as the bridge length b increases. Since 
an analytical diagonalization of A for arbitrary b is not 
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FIG. 6: Donor, bridge, and acceptor populations Po(t), 
Psif), and Pa(£) 5 respectively (from top to bottom), for 
S — 1, e/A — 5, a — 1/4, u> c /A = 5, and (top to bottom 
for bridge and acceptor, bottom to top for donor) h/3A = 
0.05, 0.1, 0.2, 0.3, 0.4, and 0.5. Error bars correspond to one 
standard deviation in each direction. Dotted lines denote the 
respective equilibrium populations from Fig. [5] of the bridge 
state. The solid lines are guides for the eye only. 
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FIG. 7: Sequential transfer rates Fob obtained from the 
bridge state (open diamonds) and donor and acceptor state 
(closed diamonds) for S = 1, e/A = 5, a = 1/4, and lo c / A = 5 
as a function of inverse temperature. The solid line refers to 
the golden rule rate 11221 . the dashed one to the Marcus rate 
Im- 



possible, we invoke an effective DBA model to capture 
the key features and compare with the numerical results. 
Accordingly, we follow Ref. [9j and consider in addition 
to the donor and the acceptor population only the bridge 
populations Pbd and Pba coming from the donor and ac- 
ceptor, respectively. The dynamics of these populations 
V = {Pd, Pbd, Pba, Pa) is then determined by V = AV 
with the reduced rate matrix 



A = 



—Fob 


rv 


r,, 







-r r - r d 














-rv - r d 







R, 




— ^DB 



(40) 



Here, T r denotes the recrossing rate from the bridge back 
to donor/acceptor and the rate for the diffusive mo- 
tion from one end to the other end of the bridge. The 
smallest eigenvalue can now be found analytically and 
reads 



r» = 



DB 



1 + fl r + f-ld - \/(5 + Mr + Md) 2 - 



(41) 

where for convenience we introduced /Lt r /d — T r /d/T db- 
In Ref. it was found that fi r /fid = b. Since Tob is 
independent of the bridge length, the ratio /id can be 
inferred from the mean first passage time of a particle 
diffusing across a one-dimensional wire [3^ |. This scales 
like b 2 so that /i^ = l/(vb 2 ) where v is a constant. Based 
on these results two scenarios for Tn(b) can be distin- 
guished. For short and moderate long bridges with b 
sufficiently smaller than 1/v (fi d > 1), Eq. (|40|l gives rise 
to r(6) « T DB 2/(1 + 6), while for longer bridges b > l/v 
(fid "C 1) the approximation Tp\b) » Tbb %/{vb 2 ) ap- 
plies. In the first case the bottleneck of the transfer pro- 
cess is the donor-bridge activation encoded in Tdb, in 




FIG. 8: Electronic populations for e/A = 2.5, a — 1/4, 
u c /A = 5, h/3A = 0.1, and 1 < S < 3 1 / 2 (6 = 1,..., 10). P { * ] 
denotes the acceptor population for a system with b bridge 
states. 
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FIG. 9: Relaxation rate Y R for e/A = 2.5, a = 1/4, 
lu c /A = 5, h/3A = 0.1, and 1 < S < 3^2 (1 < b < 10) 
(open diamonds). The solid line resembles Eq. (I42|l . The 
rates in the inset correspond to the golden rule rate 1221 . 



the latter one the rate is limited by the diffusive motion 
along the entire bridge described by 1^. An effective for- 
mula comprising these limiting cases is given by 



r B (6) 



DB 



1 



■b 2 



(42) 



Indeed, the numerical data shown in Fig. are quan- 
titatively captured by the result l|42l) with v — 0.159. 
Due to this excellent agreement between sequential trans- 
fer model and the PIMC data on the one hand and be- 
tween the corresponding relaxation rates and the expres- 
sion l|40|) on the other hand, we extracted the relaxation 
rate from the full rate matrix A for bridges with up to 
b = 100 and confirmed the validity of Eq. lj4TJ|) . 

A functional behavior of the form Tn(b) oc 1/b has been 
noted empirically already in Ref . 8] , derived theoretically 
in Ref. |2| (in a slightly different model though), and 
observed experimentally in a variety of molecular sys- 
tems I n lead- molecule-lead junctions it gives 
rise to Ohm's law in the current voltage characteristics. 
Our above findings indicate that, at least in symmetric 
DBA systems, this behavior is characteristic for situa- 
tions where T^b is sufficiently smaller than the bridge 
diffusion rate. However, a changeover to T(b) oc 1/b 2 
should be observable for very long bridges. 



VI. CONCLUSIONS 

Electron transfer processes across molecular chains 
have gained much interest in recent years, particularly in 
the context of molecular electronics. Motivated by previ- 
ous studies, in this paper the electronic dynamics across 
symmetric donor-bridge-acceptor systems has been an- 
alyzed based on a numerically exact treatment of the 
corresponding dissipative quantum system within a real- 
time MC scheme. 



Upon improving existing filter techniques we have been 
able to achieve accurate data for the population dynam- 
ics even for times where equilibration sets in or has al- 
ready been established, and for an increasing number of 
electronic bridge states. This opens the door to consider 
a variety of electronic site topologies as well as the im- 
pact of external time dependent fields. Eventually, by a 
proper reduction of the Hilbert space also the reduced dy- 
namics of continuous degrees of freedom with a discrete 
energy spectrum of the bare system are now accessible. 
Corresponding work in this direction is in progress. 

Physically, the main issue has been the role of superex- 
change and sequential transport in symmetric DBA com- 
pounds. Here, results based on Redfield theory and sim- 
plified system-bath couplings indicated different regimes 
for the length dependence of net transfer across the 
bridge comprising an exponential and an algebraic fall- 
off, respectively. The numerical MC data have been an- 
alyzed by exploiting the fact that (i) the ET dynamics 
can in the range of parameters considered be mapped 
onto Master equations with time independent transition 
rates and (ii) the path integral formulation reproduces 
the known analytical findings for these rates in certain 
limits. For the DBA compound the parameters were cho- 
sen to guarantee strong quantum effects, i.e. moderate 
dissipation, sufficiently low temperatures, and a broad 
distribution of bath modes. 

The conclusion of the numerical results together with 
the analytical findings is that in a symmetric DBA sys- 
tem (degenerate donor and acceptor energy) even for 
weaker dissipation superexchange can be expected to 
dominate only for classical or nearly classical bath. For 
quantum mechanical baths as they have been studied 
here, it does not play a significant role for ET pro- 
cesses that can be captured by time independent trans- 
port rates. Particularly, in contrast to earlier studies 
we found no dominant role of the superexchange mech- 
anism neither for relatively high lying bridges, nor for 
lower temperatures or more bridge sites. For quantum 
baths it requires extremely high lying bridge states with 
essentially frozen dynamics or temperatures that are so 
low that coherences give rise to oscillating population 
dynamics associated with a breakdown of a conventional 
rate description. Experimentally, a clear observation of a 
changeover from tunneling to hopping dominated trans- 
fer in molecular chains has not been found yet, a fact, 
that may to some extent be attributed to the scenario 
we have revealed. The results reported in Ref. show 
some indications but not a direct proof though, since the 
energetic landscape of the molecular structures changes 
with their length. By varying the number of bridge sites 
an algebraic decrease of the relaxation rate could be con- 
firmed for shorter as well as longer bridges. Based on a 
reduced DBA model we find a changeover from a range 
where the donor-bridge activation limits the transfer to a 
domain where the nondirectional diffusive motion along 
the wire is decisive. Previous simulations on biased three- 
state systems suggest in accordance with the analytical 



13 



results that an energy gap between donor and acceptor 
(in addition to the energy gap between the bridge and 
donor/acceptor) stimulates the occurrence of superex- 
change. 



alternative is given by 



%c(Ui},i7, + i) = |i^ (1) (0,6)---fc (9) fe)| (A5) 
with the real-valued matrices 
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APPENDIX A: COMPUTATIONAL DETAILS 

For the sake of computational savings, we calculate 

M 

fc = l,...,g+l (Al) 

rather than Eq. 1)26(1 . such that one single MC trajectory 
can be used to compute P s ,,si for all times tfe = (fe— 1)t, 
where the simultaneous measurement on the forward and 
backward path enhances the statistics. Equation l|38(l 
then becomes 



exp 




+ ^ CfeAfe-j 

x^>/)fe-,...,^---fcW(^,r7 9+1 ) 



(A2) 



where K^' Sf '(^k, ■ ■ ■ , £g) is defined by 



xK^ v . +1 (^,...,^) (A3) 



1 



and again is a row vector due to 771 = 2sj, while 

& q \€q,T)q+i) denotes the rjg+i's column of K^(^ q ). 

Another significant reduction of the computational 
costs of the MC calculations can be obtained by opti- 
mizing the MC weight. While 

WMc({£j},Vq+i) = expl - & A k-j£j\ 

\ k>j=2 ) 

x ^«(6,...,^)...fcW(£ 9 ,77 g+ i)| (A4) 

generally represents a natural choice for the MC weight, 
for systems with sufficiently weak damping an attractive 



KAt£) = |((r? + 0/2|ex P Hr-ff dLS /?i)|(r7' + O/2}| , 

(A6) 

which essentially neglect all dissipative effects. Although 
MC runs utilizing wmg exhibit somewhat poorer statis- 
tics than those using u>mCj the switching from complex 
to real- valued matrices which furthermore, unlike the 
can be calculated and stored before executing the 
MC moves results in a significant speedup of the simula- 
tions. In fact, for the parameters investigated here, using 
Eq. I|A5() over l|A4J> could reduce the computational costs 
by a factor of roughly four. Note that the restriction 
of the use of wmc to the weakly damped regime does 
not really impose a severe limitation since for sufficiently 
strong dissipation the overall sign problem will be weak 
enough to be tackled by any conventional PIMC scheme. 



APPENDIX B: POPULATION DYNAMICS FOR 

S = l 

For S = 1, explicit and transparent expressions for the 
electronic populations can be gained from Eq. I|18|) and 
(|19|l . Allowing some transient motion for times t < itrans, 
one easily obtains 



Pd (^trans + t) — 

trans JJ^- 



- [Pg-P B (ttrans)] 

— Pb (^trans) l , 
Pb Straus + t) = 

[Pb ~ -PbCWxs)] 



1 — exp 



(r DB +r s )t 



p, 



B 



1 — exp 



Tub 

P'OO 
R 



+ Pb (Straus) j 
Parana +*) = 
1 ' 



1 + [-Pa Straus) - Pd (ttrans)]e 



-[PS'-Psitta 
— Pb (ttrans) / j 



1 — exp 



-(r D s+r s )t 



DB 



(Bl) 
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